# maximum entropy R script that gets called by Stata
setwd("~/Dropbox/Benton-Jordan-Philips/final-JOP-replication")
library(haven)
library(meboot)

# Parse argument from rscript call in Stata
args = commandArgs(trailingOnly = "TRUE")
if (length(args)) {
    numberofvar <- args[1] # number of var to meboot
    dataset <- args[2] # dataset to call in
    T_len <- args[3] # length of T
    reps <- args[4] # number of reps
}

#numberofvar <- "1"
#dataset <- "dgp_scenario1.dta"
#T_len <- "250" 
#reps <- "1000"

me_data <- read_dta(dataset)
set.seed(13098)

# function to create nicely named vars
create.bootstraps <- function(y, first.obs, last.obs, n, x, z) {
    og_x <- x[first.obs:last.obs]
    og_z <- z[first.obs:last.obs]
    me.boot.store <- meboot(y[first.obs:last.obs], reps = n)
    the.me.boots <- data.frame(me.boot.store$ensemble)
    names(the.me.boots) <- paste0("meboot_y_", 1:n)
    time <- first.obs:last.obs
    combined.boot <- data.frame(time, y[first.obs:last.obs], og_x, og_z, the.me.boots)
    combined.boot
}

# take ME draws
me_data <- create.bootstraps(eval(parse(text = paste("me_data$y_", numberofvar, sep = ""))),
                             first.obs = 1,
                             last.obs = as.numeric(T_len),
                             n = as.numeric(reps),
                             eval(parse(text = paste("me_data$x_", numberofvar, sep = ""))),
                             eval(parse(text = paste("me_data$z_", numberofvar, sep = ""))))

# write to csv to open in stata
write.csv(me_data, file = 'me-draws.csv')
